Lots of data in our field(s) are inherently spatial
R has lots of tools for interacting with and analyzing spatial data
Lots of data in our field(s) are inherently spatial
R has lots of tools for interacting with and analyzing spatial data
sp packageThe sp package is the workhorse of spatial mapping and analysis in R.
It defines spatial classes (e.g. SpatialPoints, SpatialPointsDataFrame) that make it possible to work with spatial data.
We often just have coordinate data, but this works differently in R from fully spatial data
myCoords <- data.frame(long = runif(20, min=35, max=36), lat = runif(20, min = 3, max=5)) plot(myCoords)
You can use the SpatialPoints() function to make your coordinates into a fully spatial object
Warning: This function will assume that your columns are in the order: X coordinate, Y Coordinate
For a gold star: Which is the x and which is the y when we are dealing with latitude and longitude?
library(sp) spMyCoords <- SpatialPoints(coords = myCoords)
Now there are a plethora of spatial things you can do with these points that you couldn't before
bbox(spMyCoords)
## min max ## long 35.009496 35.923433 ## lat 3.079992 4.984301
plot(spMyCoords, axes=TRUE)
The maps package provides basic maps that can serve as a backdrop to your points
library(maps) world <- map(database = "world")
USA <- map(database="state", fill=TRUE, col=c("red","blue"))
kenya <- map(database = "world", region="Kenya")
Plot the spCoords points as solid blue filled circles on top of the Kenya map
library("RgoogleMaps")
UNCG <- GetMap(center = c(36.0674908,-79.8098429), zoom = 17)
PlotOnStaticMap(UNCG)
UNCGSatellite <- GetMap(center = c(36.0674908, -79.8098429), zoom = 17, maptype = "satellite") PlotOnStaticMap(UNCGSatellite, lon=-77.0493, lat=38.899, pch=16, col="red", cex=2)
Get a satellite image of Africa, and plot a blue + symbol on Khartoum
You can do almost anything in R that you can do in ArcGIS.
This is accomplished through additional packages, notably rgdal which is an R interface to the Geospatial Data Abstraction Library (GDAL), which is a great piece of open source software for spatial analysis. GDAL is the brains behind QGIS.
The raster package deals with geospatial raster data, relying heavily on GDAL, as well.
The spatstats package implements a variety of spatial statitics.
rgdalrgdal can read in most any vector data
The readOGR() function reads in vector data
First, download the shapefile from http://uncg.rstats.ninja/datasets/spatial/africa_simple.shp.zip. Then uncompress the file.
library(rgdal)
africa <- readOGR("/Users/wabarr/Dropbox/UNCG-workshop/datasets/spatial/africa_simple.shp",
layer="africa_simple")
## OGR data source with driver: ESRI Shapefile ## Source: "/Users/wabarr/Dropbox/UNCG-workshop/datasets/spatial/africa_simple.shp", layer: "africa_simple" ## with 1 features ## It has 3 fields
rgdalspatial data can be plotted with the base R plot() function
plot(africa, col='maroon')
raster packageFirst, download the raster from http://uncg.rstats.ninja/datasets/spatial/africa_landcover.tif
library(raster)
landCover <- raster("/Users/wabarr/Dropbox/UNCG-workshop/datasets/spatial/africa_landcover.tif")
plot(landCover)
proj4string(landCover)
## Loading required package: raster
## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
proj4string(africa)
## [1] "+proj=eck6 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
Notice that the africa shapefile is in a different projection (Eckert VI equal area) rather than latitude / longitude, which is the most common coordinate system you will see
You can reproject vector data with sp::spTransform()
You can reproject raster data with raster::projectRaster()
Let's reproject the africa shapefile into lat/lon coordinates
africaLatLon <- spTransform(africa, proj4string(landCover))
Now that they are in the same coordinate system we can plot them together
plot(landCover) plot(africaLatLon, add=T, border="purple", lwd=5)
maskedLandCover <- raster::mask(landCover, africaLatLon) plot(maskedLandCover)
sites <- read.table("http://uncg.rstats.ninja/datasets/spatial/africanSites.txt",
header=T, sep="\t")
sitesSP <- SpatialPointsDataFrame(coords=sites[,1:2], data=sites["siteID"])
plot(maskedLandCover)
plot(sitesSP, add=T, pch=16, col="blue")
Make 1 degree buffers using the `rgeos:gBuffer
rgeos::gBuffer(sitesSP, byid = T, width=1)
## class : SpatialPolygonsDataFrame ## features : 5 ## extent : -5.509803, 27.91958, -22.9821, 27.97751 (xmin, xmax, ymin, ymax) ## coord. ref. : NA ## variables : 1 ## names : siteID ## min values : site1 ## max values : site5